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Abstract We present a stochastic, time-discrete boolean model which mimics the meso- 
scopic dynamics of the desorption reactions A + A — > A + 5* and A + A ^ S -\- S in a. ID 
lattice. In the continuous-time limit, we derive a hierarchy of dynamical equations for 
the subset of moments involving contiguous lattice sites. The solution of the hierarchy 
allows to compute the exact dynamics of the mean coverage for both microscopic and 
coarse-grained initial conditions, which turn out to be different from the mean field pre- 
dictions. The evolution equations for the mean coverage and the second order moments 
are shown to be equivalent to those provided by a time-continuous Master equation. The 
important role of higher order fluctuations is brought out by the failure of a truncation 
scheme retaining only two-particle fluctuation correlations. 

PACS number (s): 05.40.-a, 05.45.-a 



I. INTRODUCTION 

When reactive processes take place in low dimensional spaces, the mean field (MF) 
laws of classical kinetics, in which each particle is assumed to interact with the system 
as a whole, become questionable. In this simple one-particle picture, information about 
special geometric constraints or many-particle correlation effects is absent. However, 
these features become increasingly important as the dimensionality of the space decreases 
and may eventually give rise to a crossover between universal macroscopic MF behavior 
and specific, strongly lattice-dependent behavior below a critical dimension d c . The 
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critical dimension d c is found to depend on the characteristics of the support (such as 
average coordination number) and the stoichiometry of the reactive dynamics, i.e., its 
degree of nonlinearity. 

The crossover phenomenon mentioned above is not specific of reactive systems, but 
is common to many other statistical systems (interacting spins, random walks, etc.). 
To understand its real nature, it is necessary to resort to a more refined description in 
which relevant fluctuation effects in concentration and occupation number space can be 
accounted for. In view of the difficulties encountered by a full-scale microscopic analysis, 
one typically tries to devise "minimal" models which contain the necessary ingredients 
to observe a significant departure of MF behavior at sufficiently low dimensions. A 
probabilistic, time-continuous approach of this kind for a detailed study of fluctuations is 
provided by Master equation (ME) techniques, in which the local dynamics is specified by 
transition rates between different microscopic states. In the ME approach, one is typically 
led to a set of evolution equations for the moments of the probability density, which can 
be taken as a starting point to compute the most characteristic macroscopic observables, 
generally associated with the first low-order moments. However, if the transition rates 
are nonlinear, the low-order moments will not obey closed equations but will rather be 
linked to the higher order ones by an infinite hierarchy of coupled equations. In the 
general case, these hierarchies cannot be solved in a closed form. Suitable truncation 
techniques may then help to overcome this difficulty. 

A more recent, very successful method to study the evolution of low dimensional 
reactive systems is based on Monte Carlo (MC) simulations [|J, in which the microscopic 
state of the system is updated at discrete time steps according to a given dynamical rule. 
One of our goals in this paper is to stablish a link between the physics underlying this 
simulational approach and the ME formalism. To this end, we shall formulate mathe- 
matically the evolution law prescribed by the MC algorithm as a stochastic dynamical 
system whose state variables are Boolean occupation numbers for the state of each lattice 
site @||]. As we shall see, this cellular automaton (CA) model yields a set of moment 
equations similar to that obtained from a time-continuous ME. 

So far, most of the literature devoted to the effects of dimensionality on reactive 
dynamics has focused on diffusion-controlled reactions j§-|| . In these systems, deviations 
of MF behavior at low dimensions may be expected due to the reduced effective mobility 
of the reactants. There are, however, situations in which diffusion can be neglected 
within the time scale of interest, like e.g. chemical reactions in solid materials and 
certain radical isolation problems . Such systems of immobile reactants constitute the 
object of our study in this paper. Clearly, the absence of diffusion reduces the ability 
of each particle to interact with all the others, thus driving the system further away 
from the applicability conditions of MF theory. This effect is enhanced in the presence 
of short range interactions restricted to, say, nearest neighbours and hard core exclusion 
not allowing more than one particle per lattice site IHJ. 

More specifically, we will study the ID lattice dynamics in two particular examples 
of nonlinear irreversible associated with the cooperative desorption systems (CDs) 

A + A A + S cooperative partial desorption (CPD) (la) 

A + A — ^ S + S. cooperative total desorption (CTD). (lb) 

where A is the reactive species, S the empty lattice site and is the rate of reaction. 
Notice that steps (pTa )- (|T5 ) are typical parts of realistic, more complex reaction schemes 



as discussed further in the Conclusions. Despite their apparent simplicity, they display 
a complex non-MF behavior characterized by frozen, non universal, initial-condition- 
dependent steady states [0,0,0]. The dependence on the initial conditions is the 
signature of the weak ergodic properties of the CDs, which is in turn related to the 
irreversibility of the reactive schemes flIaD-(Ib"). A remarkable feature of the CDs is 



the special structure of the underlying moment equations, which allows one to compute 
explicitly the coverage and the fluctuation dynamics and thus test in a very efficient way 
the applicability of diverse truncation approaches. 

The paper is organized as follows. In section [TT], a comprehensive study of the one- 
dimensional CPD is presented. First, we introduce a particular biased implementation 
of the CPD and characterize its non-MF behavior using an intuitive argument. In [II B| , 
we introduce a stochastic dynamical rule which mimics the microscopic, time-discrete 
evolution of the system as prescribed by MC simulations and show that this dynamical 
rule characterizes properly the steady states of the system. In [11 C] , we derive evolution 
equations for the first- and second-order moments of the probability distribution for the 
microscopic states and use the same formalism to obtain a hierarchy of coupled equations 
describing the time evolution of a particular subset of moments, namely those involving 
contiguous occupation numbers. These moments can be identified with the probabilities 
of finding a randomly chosen string of contiguous sites simultaneously occupied. In |II D| , 
it is found that the solution of the hierarchy for both microscopic and coarse-grained 
initial conditions leads to an explicit expression for the mean coverage of the lattice. 
We also derive equations for the coarse-grained dynamics of vacant sites in terms of 
generalized probabilities for occupied sites. Next, in |11E| , we present a general ME 
for an array of Boolean variables and show that, in the special case of the CPD, the 
evolution equations for the mean occupation number and the second-order moment are 
similar to those obtained in the CA approach. In section |HI| , we study the relevance of 
inhomogeneous fluctuations in the CPD system. A set of evolution equations retaining 
only pair fluctuation correlations fails to reproduce the correct dynamics, the reason being 
the nonnegligibility of higher order fluctuation correlations. These can be calculated 
explicitly from sets of generalized moment equations, which also allow us to compute the 
dynamics of the vacant sites for clusters of small size. In section [IV], we extend the main 
results of sections [TI| and |TTT| to the CTD. The conclusions are summarized in section [V], 
which also outlines a research strategy for future work in this area. In the Appendices, 
we deal with a symmetric version of the CPD and discuss its equivalence with the biased 
CPD. We also study the dynamics of particle islands. 



II. BOOLEAN DYNAMICS OF THE CPD: MODEL AND EXACT RESULTS 

A. Definition of the system and MF rate equation 

Consider a ID lattice with N sites which may be either empty (S) or filled with a 
single particle (A). The reaction takes place according to the scheme of Fig. [I]. At 
each time step At, a lattice site is randomly chosen. If the selected site is filled, the 
particle desorbs with a certain probability (equal to the rate of reaction kn) provided 
that its right neighbour site is also occupied; otherwise, nothing occurs. The classical 
rate equation for this CPD scheme (Eq. (0)) reads 



(2) 



where c(t) is the A particle concentration, here identified as the mean fraction of occupied 
sites, also referred as "coverage" . Its solution predicts a decay of the form 

to a zero-coverage steady state. However, as we see below, the restricted ID geometry 
of the lattice allows for non MF steady states. It is indeed easily seen that all lattice 
configurations without contiguous particles cannot evolve in time. To understand this 
obvious MF failure in detail, we need to set up a model for the microscopic dynamics. 



B. Boolean modeling and steady states 

We stipulate that each site in the lattice is characterized by a Boolean occupation 
number: 

f 1 if site i is occupied . , T , 

Tk(?t) = { n ...... : , 1 = 1,. ..,n. (4) 

[ (J if site i is empty 
Each site evolves according to the nonlinear dynamical rule 

n i (t + At) = n i (t)-$(t)Z R (t)n i (t)n i+1 (t), i = l,...,N (5) 



where 



( &\t) \ 



$(f) = 0,l. (6) 



is a N- dimensional vector of stochastic Boolean decision variables. If, say, site % is chosen, 
the 2-th component will be one and all other components will be zero. Hence, the different 
components £jy (t) in a given, single realization are correlated: 

N 

£d?(*) = l (7) 

8=1 



Since at each time step At, one site out of N is randomly chosen, the mean value £jv (£) 
over an ensemble of realizations will be 1/N. The additional stochastic variable £r(£) 
takes randomly the values 1 and with probabilities p R and 1 — p R , respectively. 

To have a well posed problem, we still have to specify the boundary conditions. We 
shall take periodic boundary conditions, which amounts to setting n N+ i(t) = ni(t) in 
the last equation (|5|). 

Implicit in the above scheme is the idea that, at most, one reactive event may occur 
in the lattice in a given time interval At, even if there are several a priori reactive 
pairs. The justification of such a prescription is that, typically, the ocurrence of a 
reaction requires to overcome an activation energy threshold. At ordinary temperatures 



(supposed to prevail here), this can be achieved only if a sufficiently strong fluctuation 
impinges on the system as a result, say, of the coupling between the adsorbate and the 
lattice. Since the probability of such an event is small, one may expect that reactions 
will occur asynchronously at different lattice sites, which is precisely what is stipulated 
in our evolution rule. 

Eqs. (H) will be used in the sequel to describe the dynamics of the CPD at three 
different levels: 

1. A "microscopic" level in which the exact microscopic dynamics, i.e. the ini- 
tial condition {^-i(O)} and the decision path represented by the variable sets 
{£/v(0), . . . , £tv(£)} and {O?(0), • • • £,R,(t)}, is assumed to be known in detail. 

2. an intermediate level in which the system is prepared in such a way, that the initial 
condition is known in detail and one averages only over the subsequent dynamics, 
i.e. , over the different paths or realizations. 

3. a "coarse-grained" description in which one averages over a statistical ensemble of 
realizations and a nonequilibrium ensemble of initial conditions. 

After a time t st a steady state {nf } is attained. The dynamical rule (|5|) allows a straight- 
forward characterization of these steady states. Since the first term in the r.h.s. of (|5D 
cancels with the l.h.s., one has: 

d?W«i = 0. (8) 

This holds for all times t > t st . Since £jv (0 is a random variable, it will be nonzero at 
some time t for any given site i, meaning that the product nf nf +1 must vanish for all i. 
This expresses the fact that two contiguous sites cannot be simultaneously occupied in a 
steady state. However, not all configurations satisfying this condition are attained with 
equal probability. For instance, if N is even, we can be certain that the steady-state 
configuration with alternating occupied and empty sites will never be reached if we start 
with less than half the lattice filled. Therefore, a simple combinatorial counting of the 



non-evolving configurations fails to provide the correct mean steady-state coverage WA 



C. Moment equations and cluster dynamics 

Let us now turn to the dynamics generated by Eqs. @. First, we fix the initial 
configuration (ni(0), . . . ,n N (0)) of the lattice and take the average of Eq. (|5|) over an 
ensemble of different realizations. Using the statistical independence of 6v(£), £r(£) and 
the occupation numbers ni(t), this yields 

Vr 

n i (t + At)-n i (t) = -^n i {t)n l+1 {t), i = l,...,N. (9) 

The quantity ^3 can be regarded as the probability that a reaction takes place at a given 
site % in the time interval (t, t + At) given that site i and i + 1 are occupied. The rate of 
reaction is obtained by dividing this probability by the time step At: 



' R NAt 



(10) 



To derive the continuous-time limit of Eq. (|9]), we divide both sides by At and let this 
quantity go to zero in such a way that the reaction rate k R remains finite. This can 
be e.g. accomplished by letting the system size N simultaneously go to infinity so that 
NAt = C . We shall choose the constant C equal to unity, implying that, after one time 
unit, each lattice site has been visited once on average. With this choice, k R becomes 
numerically equal to the value of p R . The resulting equation reads 

d 



—n i (t) = -k R n i (t)n i+l (t) (11) 
at 

In a similar way, one can obtain an evolution equation for the second-order moment 
ni(t)rij(t). Multiplying (|5|) by rij(t + At) and using (^), one finds 

PR —7T\ 7T\ 7TT PR 



ni(t+£t)nj(t+At)-ni(t)nj(t) = -— ni(t)n i+l (t)nj(t) - —m(t)nj(t)nj +1 (t), i ^ j 

(12) 

which in the continuous-time limit becomes 
d 



ni(t)rij(t) = -k R ni(t)n i+1 (t)nj(t) - fc R n i (t)n i (t)n J+1 (t). (13) 



For j — % + 1, Eq. ( p~3"D together with the Boolean property nf{t) = rii(t) yields 

— ni(t)n i+1 (t) = -k R rii(t)n i+1 (t) - k R ni(t)n i+1 (t)n i+2 (t) (14) 

Unless otherwise specified, we shall set p R = 1 for simplicity, implying that k R is equal to 
unity. Note, however, that it is possible to absorb the reaction rate k R appearing in the 
time-continuous equations (pTT|) and (|l3f) in the time scale by means of a dimensionless 
variable r = k R t |T0|JT4]1 . 



An important feature of Eqs. (|Tl|) and (|T4{ ) is that they couple in a linear fashion 
moments involving only contiguous sites. This turns out to be a generic property that 
also holds for the dynamics of higher order moments. 

Let us now consider a more general subset of k contiguous lattice sites i, . . . , i + k — 1. 
The probability to find all sites occupied at time t, i.e., a cluster of size k, is given by 
the quantity 



Mj?{t) = 1 n \j{t\ i, k = 1, . . . , N. (15) 

j=i 

To derive an evolution equation for Mj?\t), we must generalize Eqs. (0) and (|I~4"D. 
Taking the dynamical rule @ as a starting point, one obtains 

i+k-l i+k-l i+k-2 i+k-1 i+k 

n nj {t+At)-i[ ni(t) = -E $(*)£«(*) n ^w-^^'w^wn^w. 

j=i j=i j=i j=i j=i 

(16) 



where % — 1, . . . , N. Taking averages in ([If]) for k < N, one gets 



Mf (t + At) - M®(t) = -^j^ Mf (t) - 1 MtUt). (17) 
We now introduce the global quantity 

i JV 

H N) (t) = j,J2M^(t), (18) 

representing the probability of finding k filled adjacent sites in the iV-site lattice. Sum- 
ming over the site index i in Eqs. (|ITD , we find 

pW( i + A i )-Pr (i ) = -^pf , W-^SW, k<N. (19) 

Taking the continuous-time limit of (|1^) for a given finite value of k, a set of evolution 
equations for Pk{t) = linijv^oo Pk N \t) is obtained: 

dP (t) 

^i = _(A;_i) j p fe (t)- j p fc+1 (t) 5 A; = 1, 2, ... , k max , (20) 

where the integer k ma x is the size of the largest cluster. The first term on the r.h.s. of 
(p0|) corresponds to a destruction of a k-particle cluster by the interaction of two particles 
inside the cluster, while the second term stands for its destruction by desorption of its 
rightmost particle, which requires an additional occupied site, i.e., a cluster of k + 1 sites. 

If the lattice contains only clusters of finite size, then k max < oo and the hierarchy 
of equations (|2H| ) is truncated by the condition Pk max +i(t) = 0. The last equation ( |2T)| ) 
reads then 

dPk J {t) =-(k max -l)P krna M- (21) 



D. Solution of hierarchy and mean coverage for microscopic and coarse-grained 

initial conditions 

Suppose that we start with a given initial configuration {^j(0)} characterized by the 
nonvanishing set of probabilities 

{A(o),A(o),...,/\_(o)}. (22) 

Using the linearity of Eqs. (|20|), one can show that the solution for t > is of the form 

k ma x-k C exD ('_^ > ) —1)3 

P k (t)=exp(-(k-l)t) £ 1 Pl . -n+i(0), k = l,...,k max . (23) 

3=0 I' 

Note that for t — > oo only Pi survives, representing the probability that a randomly 
chosen site be occupied. Clearly, this probability is equal to the ratio between the 
average number of occupied sites and the total number of sites, i.e., the mean coverage 
of the lattice c(t). In the long-time limit, one has: 



c(oo) = ^(oo) = V^^pl= £ [ -f-P j+ i(0). (24) 

We see that the asymptotic value of the coverage depends on all details of the initial 
cluster distribution. The validity of Eq. ( [2"4] ) is confirmed by MC simulations (Fig. Q(a) 
). The simulations were performed according to the time-discrete rule @. The data 
correspond to a chain ring of 100 particles and the initial configurations: 

CI = 11111111111111110111111111111011110111111111111111011111111111111111 

11101111111011011111111011111010 
C2 = 01111001111110010011110111110111011111111110100011111111011011001111 

11011101001111101011111111101111 
C3 = 11010101101001100100101010010010000101001011101001001001011001001111 

01100011111001011001101001100111 
CA = 10011100111000010011011000001000000000010000100100100100011000101000 

01001010010101000011001101100010 

Let us now assume the more realistic situation in which the initial configuration of the 
system is not known in detail, but is characterized only through the particle coverage p. 
If we consider a uniform ensemble of all random microscopic configurations with a given 
p, the system will initially be translationally invariant and this property will propagate 
in time. Thus, the probabilities of finding k consecutive filled sites i, . . . , % + k — 1 

Pl N \t) = ( + f[ 1 n J (t)}, k = l } ...N (25) 
j=i 

(the angular brackets (• ■ ■) stand for the averaging over realizations and initial states) will 
not depend on the evaluation site i. After averaging Eq. (|T7|) over the initial conditions 
and going over to the time continuum, one obtains: 

dP (t) 

-J±l = —(k — 1) P k (t) - P k+1 (t) k = 1, 2, ... , (26) 

which is formally identical to Eq. (|20|) for P& with k max — > 00. However, the initial 
conditions are now simpler. Since the initial occupancies of each site are statistically 
independent, one has: 

P k (0)=p k , A; = 1,2,... (27) 



A solution satisfying the initial conditions (|27j ) can be easily found by using an ansatz 
suggested in []10| . One obtains 

P 1 (t) = p exp [p (exp (-t) - 1) ], (28a) 
Pk(t) =p k - 1 ex P [-(k-l)t]P 1 (t), A; = 2,3,... (28b) 

Again, the survival expectancy Pi(t) corresponds to the mean coarse-grained coverage 



c(t) = limjv^oo ■ F° r short times, the series expansion of (|28a|) yields: 



c(t) =P 1 (t)=p-pH + /(P 2 + 1} t 2 + o(t 3 ) (29) 

It is interesting to compare this with the short-time asymptotics of the exact solution of 
the MF equation @ for &r = 1: 

CMF(t) = -^-=p-p 2 t + pH 2 + o(t 3 ) (30) 
l+pt 

Thus, we see that the first MF deviation appears in the second order term. For long 
times, one observes an exponential decay into an absorbing state characterized by the 
asymptotic coverage 

c(oo) = p exp (— p) (31) 

Formulae (|28|) for the transient behavior and (|31|) for the asymptotics are confirmed 
by time-discrete MC simulations (Fig. |2](b) ). As in the previous case, all Pk{t) with 
k > 2 vanish for t —>■ oo. Note that Eqs. ( p8|) are formally obtained from Eqs. (^) by 
averaging over a uniform distribution of microscopic states with global coverage p. This 
amounts to allowing for clusters of arbitrarily large size , i.e., k max — > oo. 

So far, we have only considered the dynamics of particles, but it is also of interest to 
ask for the time evolution of empty intervals of a given size inside the lattice. Clearly, 
the number of empty intervals must grow as the reactions proceed. To set up explicit 
equations for clusters of vacant sites ("holes"), we consider the time-discrete evolution 
equation for the complementary occupation number Si(t) = 1 — rii(t): 

s i (t + At)-s i (t)=^ ) n l (t)n i+1 (t), i = l,...,N. (32) 

By averaging this equation over an ensemble of realizations and a uniform initial distri- 
bution of microscopic states with a coverage p one gets: 

^ = W). (33a) 

where Si(t) = (si(t)}. Eq. ( |33a| ) is in agreement with the conservation law + 
Si(t)) = 0. 

The dynamics of hole multimers can be derived starting from the evolution equation 

i+k—l 

for the product J][ s&(£). The general form of the evolution equation for Skit) = 

k=i 

i+k—l 

( TJ Skit)) becomes more and more complicated with increasing cluster size. The next 

k=i 

three equations read 

^ = P 2 (t)-P 3 (t), (33b) 
^ = P 2 (t) - P 3 (t) + P 4 (t) - Pg(f), (33c) 

' /s ' ,m = p 2 (t)-Ps(t)+P4(t)-p 5 (t)-p^(t)+pm)+p^(t)-pmt), (ssd) 



dt 



where the quantities P^ m (t) on the r.h.s. of Eqs. (|33cj )- (|33d|) are the joint probabilities 
of finding two fully occupied strings of length n and m separated by I lattice spacings. 
The initial conditions for which these equations must be solved are 

S k (0) = (l-p) k , k — 1,2,... (34) 

Eqs. ( |33a| ) and (|33b|) can be integrated straightforwardly, since the functions on the 
r.h.s. are already explicitly known. In section [ITJ, we show how to obtain explicit 
expressions for the probabilities Ss(t) and S^(t) by solving a generalized hierarchy for 
the joint probabilities P^' m {t). 



E. Comparison with the ME approach 

An alternative way of deriving moment equations in the framework of a time- 
continuous approach is provided by the Master equation (ME). Let us again consider 
a set of N Boolean variables {n} = (m, . . . ,Un) arranged on the sites of a ID lattice 
with periodic boundaries (n N+1 — n\). Starting from an arbitrary state, the evolution 
of the probability density P({n}; t) is governed by a ME of the form: 

=-!>({"} "> M.Wit) + E« - {n},t)P({n'};t) (35) 

where {ri} = (ni, . . . , 1 — rij, . . . , njv) and Wj({n} — > {n'}, t) is the transition rate from 
the state {n} to the state {ri} at time t. To obtain an evolution equation for the average 
occupation number 

(n l (t)) = Y,n l (t)P({n},t), (36) 

W 

we multiply Eq. (ft5Q by rij(t) and sum over all different configurations {n}. This yields: 

j t (m{t)) = ( [1 - 2m{t)] w t ({n} - {n'}, t) ) (37) 
For the time evolution of the second order moment one gets 

- {n,(t)n 3 (t)) = ( [ nj (t) - 2n,(t) nj (t)] Wi ({n} ^ {ri}, t) ) 

+ ( [m(t) - 2m(t) nj (t)} Wj ({n} -> {ri}, t) ) (38) 
For the CPD, the transition rate reads 

Wj{{n} {n'},t) = k R ni(t)n i+1 (t). (39) 
Using Eqs. (37) and (|38f) and the transition probability (|39|) with kn — 1, one can for- 



mally recover Eqs. (JTT|) and (|T^) averaged over an ensemble of random initial conditions. 



As expected, this also holds for the higher order moments [[14]. Note that in the ME 
approach the double averaging over realizations and initial conditions is performed si- 
multaneously via the probability distribution P({n},t). In the CA model, this has been 
done in two separate steps. One first chooses randomly an initial condition with a given 
p and thereafter performs a series of different realizations. Next, one again selects a new 
random initial condition and a new set of realizations is carried out, and so on. 



III. TRUNCATION SCHEMES AND ROLE OF FLUCTUATIONS IN THE 

CPD 



The statistical properties of the N-site CPD system are described by the "Boltzmann- 
like" distribution functions c^j ...j n (i) = {n^it) n i2 (t), . . . , n in (t)}, which can be ex- 
pressed in terms of the mean occupations c\ (t) = (n>i(t)) and a set of correlation 
functions 

2 ,..., in (ti;...;t n ) = (Sn h (t 1 )6n i2 {t 2 )...Sn in (t n )), (40) 

where 5rii(t) = n^t) — c\ N \t). In the following, we shall omit the superscript N in the 
quantities c and / for the sake of notational simplicity For n = 2, and t\ = t 2 = t, Eq. 
(fiO|) defines the equal-time pair correlation function fij(t) = fij(t;t). 

An approximate closure to the hierarchy equations can be obtained by performing 
an expansion in terms of correlation functions and keeping them up to a certain or- 
der, neglecting the higher-order ones. In this so-called Ursell expansion one defines the 
correlation functions / through the relations 

Ci,j(t) = Ci{t) Cj {t) + f i!:j (t), (41a) 
Ci,j, k (t) = Ci (t) Cj {t) c k (t) +c k (t) fij(t) + ait) f jtk (t) + Cj (t) f itk (t) + fi, jlk (t), (41b) 

etc. 

where the quantities /jj, fi j t k account for the fluctuations: 

fij{t) = (Smtyrijit)) (42a) 
fijAt) = (8ni(t)6nj(t)6n k (t)), (42b) 
etc. 

Different approximations emerge from Eqs. (|4l| ) when terms containing correlation func- 
tions of higher orders are neglected. 

In the zeroth approximation all correlations between occupation numbers are ne- 
glected and the state of the system is completely specified by the mean occupation 
numbers Ci(t). From Eq. @ we have the MF equation 

Ci(t + At) = Ci{t) - ^a(t)c i+1 (t). (43) 

In the first approximation the rate equation ([f3|) is extended with a term linear in the 
NN pair correlation /j j +1 (t), 

Cl (t + At) = a(t) - ^Ci(t)ci +1 (t) + -^/ M+ i(t). (44) 

Eq. ( ^4[) describes the corrections to the MF equation ( ^3|) caused by the fluctuations 
fi ti+ i(t) that are built up by sequences of correlated reactions between the particles. 
These correlations are calculated from the second hierarchy equation (|1"2"D with j = i + 
Neglecting terms in o(l/N 2 ) and higher in the r.h.s., we have 

/ M+1 (t + At) = -l Ci+1 (*)/ i)i+2 (t) + f i>i+1 (t) - l( Q (t) Q+1 (t) - Ci (t)[c i+1 (t)} 2 

+ c i+2 (t) / M+1 (t) - Cj+i(t) fi, i+ x(t) + f i4+1 (t) ) (45) 



Here the first term in the brackets represents the propagation to neighboring sites of 
fluctuations correlated over a distance of two lattice spacings. The other terms essentially 
contain the correlations created or destroyed by the reaction process; the terms linear in 
fi t i + i represent the effect of existing pre-reaction correlations, and the quantity 

- ^ Ci(t) c l+l (t) + 1 a(t) [c t+1 (t)] 2 (46) 

accounts for the correlations created by the reaction from a completely factorized state. 
The other pair correlations are calculated in a similar way from (|T2|) with j — i + 1: 



f i>i+l (t + At) = f t (t) - - c,-(t) f i+hi+ i{t) - ^ (c i+1 (t) + c i+l+1 (t)) f iiW (t) 

- Ci+i(t) fi,i+l+l(t) + o(l/ N 2 ) + higher order terms, 1 = 2,..., [N}/2. (47) 

In the spatially homogeneous case, the system is translationally invariant; the solution 
obeys the relation q(£) = c(t) and the correlations fij(t) = f\i_j\(t) depend only on the 
relative distance between the evaluation sites. In the limit N — > oo, Eqs. (@), (|45|) and 
(0) yield the following set of time-continuous equations 

,/r(n ^-c(t) 2 -A(t), (48a) 



dt 

dm 

dt 
dfijt) 
dt 



c(t) 2 ( C (0-l)-/i(*)-c(0/ a (0, (48b) 
-c(t) [/,_! (*) + 2h (t) + f l+1 (t)i 1 = 2,3,... (48c) 



Eqs. ( ^8[) can be used as a starting point for a truncation based on the neglect of pair 
correlations f\(t) for all / equal or larger than a certain cutoff integer l c . For l c = 2, i.e., in 
the pair approximation, the system of two equations approaches the steady state (0,0), 
although the asymptotic decay of c(t) is much slower than the t~ l decay prescribed by 
the MF approach due to large transient effects in fi(t) (Fig. |3|). A log-log plot of the 
numerical solution suggests an inverse power law behavior of c(t) for long times with an 
exponent less than one (see inset in Fig. ^). In the pair approximation, the stability of 
the MF fixed point can be proven by geometric arguments Q Numerical results suggest 
that this property holds for any arbitrary l c , meaning that the system truncated to any 
order will always attain the MF steady state of zero concentration, in contradiction with 
the non-MF results reported above. This reflects the non-negligibility of the three-point 
fluctuation correlations fi,j,k(t) which, as we shall see presently, indeed become of the 
same order of magnitude as the pair correlations fi(t) already for small I. 

One can compute exactly these correlation fluctuations by considering the evolu- 
tion equation for the joint probabilities Pf\{t), which in our model correspond to the 
distribution functions 

i+j-l i+j+l+k-2 

( II ri r {t) II *.(*)>■ (49) 

r=i s=i+j+l— 1 



*This is illustrated by an analysis of the 2D velocity field 



Taking Eq. (|T6| ) as a starting point, a tedious but straightforward calculation shows that 
Pfl{t) evolves according to 



dPflit) 



-U - l)P$(t) - (k - V)Pf>{t) - P£#(t) - P% +1 (t), j, k,l = l, 2, 



>(0 



(i-i), 



,(0 



dt 



(50) 



The first (second) term on the r.h.s. of ( j50|) corresponds to an event in which a particle 
inside the cluster of j (k) occupied sites desorbs by interaction with a neighbour inside 
the cluster, while the third (fourth) term stands for the desorption of the rightmost 
particle inside the cluster by reaction with an occupied neighbour site at its right. Eq. 
(j5C|) is solved using the ansatz 



pfl {t) = p>+^ exp [- (j + k - 2)t] P?j (t) , j, k 
This leads to the (infinite) closed coupled set of equations 

dpfUt) 



(0, 



1,2, 



dt 



-pexp(-t) p$(t)+pi:r\t) 



>(!-!)/ 



Z = 2,3, 



(51) 



(52) 



with the boundary condition Pxi(t) = -Pj^)- An explicit expression for -Pii(t) can be 
obtained from Eq. ( p2[ ) by means of a generating function [|Tj|. The final result reads: 



p?Ut) = p 1 (t) 



l-l 

k=0 



111 



'Pi(ty 



V 



/k\ 



In 



'my 



p 



(53) 



This allows to compute the fluctuation correlation fi(t), which is given by 

fi(t)=P$(t)-[Pi(t)] 2 , 1 = 1,2,... 



(54) 



This function decreases superexponentially with increasing distance I, approaching zero 
for I — > oo. The three-point NN fluctuation correlation h(t) = fi,2(t) = ( <5n, Srii + i 5rii + 2 ) 
can be expressed as 



h(t) = Ps(t) - [P!(t)] 3 - (2/x(t) - f2(t))Pl(t)- 



(55) 



Formulae (|54D and ( pq) are confirmed by numerical simulations (Fig. ||). These also 
show that h(t) may already get larger than fi(t) for 1 = 2. In fact, the asymptotic values 
h(oo) and ^2(00) are of the same order of magnitude in the whole range of the initial 
coverage p (Fig. |5]). Higher order fluctuations become comparatively large in the dilute 
limit p< 1. 

In the light of the above, the failure of the truncation scheme seems thus to be due 
to the fact that it includes pair correlations of arbitrary long range whithout taking 
relevant higher order fluctuations into account from the very beginning. Although we 
have only shown the nonnegligibility of correlations up to the third order, some numerical 
computations suggest that, regardless of the initial conditions chosen, it is necessary to 
take into account the whole set of fluctuation correlations to obtain a deviation from the 
MF steady state. 



A more restrictive truncation scheme consists in neglecting the effect of large clusters 
in Eqs. ( ^6|) by setting P&(t) = beyond a given size [fTi}| . This truncation procedure 
yields a smooth expansion of the exact results and fl3T|) in powers of p. 

To conclude this paragraph, let us note that the explicit knowledge of the joint 
probabilities Pf\{t) can be used to integrate Eqs. ( |33c| )-( |33"cl| ). The solutions of Eqs. 
(j33[) read 

Sx(t) = l-Pi(t), (56a) 
S 2 (t) = l-2P 1 (t) + P 2 (t), (56b) 

S 3 {t) = l + (-3+p- l -p 2 ) P x {t) + 2P 2 (t) - ~P 3 (t), (56c) 

Si(t) = 1 + (-4 + 3p - 2p 2 + ip 3 ) P 1 (t) + (3-p + \p 2 ) P 2 (t) - P 3 (t) + ~P 4 (t). (56d) 

3 2 6 

They are shown in Fig. |6] for the case of an initially full lattice. 



IV. DYNAMICS OF THE CTD 

The starting point is again a ID lattice with empty (S) and occupied sites (A). The 
reaction now proceeds according to the scheme depicted in Fig. |7|. At each time step, one 
site is randomly chosen. If the chosen site and its neighbour are occupied, both particles 
will desorb with probability k R . The classical rate equation reads 

Mf) -2k R c(tf, (57) 



dt 



which yields a faster decay than (g) to a zero concentration steady state, but again 
proportional to t _1 for long times. It is easily seen that the actual ID system has the 
same steady states that the CPD system, i.e. all configurations with isolated particles. 
Once more, a deviation from the MF behavior is found. 

To set up the corresponding microscopic evolution law, we must take into account 
that this time a particle at a given site i will also desorb if its left neighbour site z — 1 is 
chosen and it is occupied. One has therefore an additional contribution in the r.h.s. of 
the dynamical rule: 

m(t + At) = ni (t)-^\t) £ R (t) m(t) n i+1 (t) - dT X) W ni_i(t) rn(t), % = 1, . . . , N. 

(58) 

Like in the CPD case, we can take Eq. ( |58|) as a starting point to derive evolution 
equations for the first- and the second-order moments: 

^-riiit) = -7ii_i(t) m(t) - fii(t) n i+1 {t) (59a) 

-j t n i{t) n i+i(t) = -ni-x(t) riiit) - rii-i(t) m(t) n i+ i(t) - ruif) n i+1 (t) n i+2 (t) (59b) 
d 



—rii(t) n i+l (t) = -rn(t) n i+1 (t) n i+l (t) - m(t) n i+ i-x(t) n i+ i(t) 
dt 

-rii(t) n i+ i(t) n i+ i + i(t) - n^t) rii(t) n i+ i(t) , 1 = 2,3,.... (59c) 



A generalization of Eqs. ( |59a| ) and ( |59b| ) for strings of k consecutive sites leads again to 
a hierarchy for the global probabilities P k {t)'- 



dP (t) 

_J±± = -{k - l)P k {t) - 2P k+1 (t), k — 1,2,..., k max . (60) 

Note that the prefactor of the second term is now 2 due to the additional reactive event 
between the leftmost particle inside the cluster and a particle at its left neighbour site. 
An equation of the form (|60|) has been obtained in the context of reactant isolation 
and random sequential adsorption (RSA) models JTEyT7|] . In these models, the deposition 
of a dimer on the lattice is equivalent to the desorption of a pair of reactants and empty 
pairs of sites available for deposition correspond to unreacted pairs of reactants (see ref. 

Ill)- 

The solutions of the hierarchy (^) for a given initial configuration read 

P k (t)=ex V (-(k-l)t) £ ^ } -P k+j (0), k = l,...,k max . (61) 

(c.f. ref. [p8|) from which the mean asymptotic coverage follows straightforwardly: 



c(oo) = AM = lim J -^p> = £ LJL P, +1 (0). (62) 

3=0 J 

This result is again in excellent agreement with MC simulations (Fig. |8|(a) ). The coarse- 
grained solution of Eqs. ( j60| ) is formally obtained by setting k max = oo and Pj{0) = p> . 
This yields 

P x {t) = p exp [2p (exp (-t) - 1) ], (63a) 
P k (t) =p k - 1 exp[-(k-l)t]P 1 (t), A; = 2,3,... (63b) 

Thus, the corresponding stationary mean coverage takes the form 

c(oo) = Pi(oo) = p exp (— 2p) (64) 

Note that, contrary to the CPD case, the />dependence of the asymptotic coverage is 
no longer monotonic. Thus, as long as p > 0.5, c(oo) increases when p is decreased 
(Fig. |(b) ). The initial situation of a fully occupied lattice (p = 1) corresponds in 
random dimer deposition to an initially empty lattice. In this particular case, a famous 
combinatorial argument by Flory |L9| predicts the value e -2 for the asymptotic mean 
fraction of empty sites (in our picture, occupied sites) at jamming, in accordance with 
the formula (|B~4|). 

For the CTD system, it is also possible to derive a truncation hierarchy by introducing 
fluctuation correlations in Eqs. (p9|). Like in the CPD case, the second order hierarchy 
yields a MF steady state, in contradiction with the exact solution of (0). Once more, 
the reason for the failure is the nonnegligibility of the three-point correlations, some of 
which can be computed by solving the equation 

^ = -u - Dp>) - (k - Dpfid) - pjKto 

-Pj1i,k(t) - Pi2+iW - P!k+l(t), j,k,I — 1,2,... (65) 



which is also found in models for random dimer deposition pD| . As shown in [11 Dj , the 
correlation functions Pj^ can be used to set up explicit equations for the dynamics of 
vacant sites. 

Finally, let us mention that it is possible to rederive the moment equations (|59|) and 
the hierarchy ( |60l) from a ME analogous to Eq. (|35|). 



V. SUMMARY AND OUTLOOK 

We have developed a one-dimensional CA model which mimics the mesoscopic dy- 
namics of two cooperative desorption reactions. An important advantage of this approach 
is that one directly sees the effect of a small modification of the MC algorithm on the 
underlying moment equations. The model allows us to derive a hierarchy of linear equa- 
tions from which the mean particle coverage can be computed for both microscopic and 
coarse-grained initial conditions. We have seen that a truncation scheme retaining only 
pair fluctuation correlations fails to provide the correct behavior of these systems, due to 
the importance of higher order fluctuations. However, an alternative truncation based 
on neglecting the effect of large clusters in the dilute limit yields a smooth expansion of 
the cluster dynamics in powers of the initial coverage. In the case of the CTD, we have 
pointed out some analogies with models for dimer deposition. 

The CDs studied here are relevant both from a theoretical and from an experimental 
point of view. Despite their simplicity, they exhibit a complex non-MF behavior char- 
acterized by nontrivial memory effects and spontaneous ordering. On the other hand, 
a wide number of physical processes on surfaces involve cooperative desorption of the 
products pl|p^] . In particular, several problems in the context of exciton dynamics in 



molecular cristals p3| , ^4| , recombination of condensed-gas radicals |25j and bond forma- 



tion in polymers PE|J7|j can be mapped into the CDs. A classical example consists of 
a long polymer chain (methyl vynil ketone) formed by immobile radical groups which 
react with nearest neighbors so as to form rings when the chain is heated |7[. Clearly, 
some groups will be isolated by the reaction dynamics and will remain unreacted. In the 
CTD model, the unreacted groups can be identified with the A particles and the rings 
correspond to pairs of vacancies left upon reaction. The quantity of interest, i.e., the 
fraction of unreacted groups, plays then the role of the lattice coverage 6. 

The effect of incorporating the backward reaction step S + A — > A + A to the CPD 
has been studied in refs. |p7|Jll| ) [12| . In this case, the mixing properties of the system are 



restored and one attains a MF steady state with equal number of empty and occupied 
sites. It would be desirable to extend the present study to the reversible CTD and also 
study the influence of an additional random particle input S — > A on the CDs. 

Another possible extension of our work consists of increasing the range of the local 
interactions by allowing, for instance, interactions with next to nearest neighbours. One 
expects that the system approaches the MF dynamics as the number of interacting neigh- 
bours increases. It is certainly worth to characterize this approach in a more quantitative 
way. 

The dynamics of the one-dimensional CDs changes significantly in the diffusion lim- 
ited case P|,|2g|-pT|] . One can account for the mobility of the particles by introducing 
additional diffusion terms in the dynamical rule. For the CPD, this has been done in ref. 
52 ] for initial conditions of the form fl27|). As expected, diffusion yields an (anomalous) 



decay into a zero concentration steady state. The results in ref. []32| have been compared 
to an off-lattice solution by ben-Avraham et al. |33|. Interestingly, the on-lattice solu- 
tion displays a slower decay of the coverage for early times due to the finite propagation 
velocity of a local concentration perturbation on the lattice. A CA approach for the 
diffusion-limited CDs has been developed by Privman [33]. In the Privman model, all 
lattice sites are synchronously updated at each time step. An extension of this model 
for the diffusionless CDs studied above is also worth carrying out. 

A natural generalization of our calculations is to study the dynamics of the CDs on 
Bethe lattices of arbitrary dimension, for which results for the cluster dynamics derived 
heuristically by Majumdar and Privman (see |H|) are available. For the CPD, we expect 
to obtain results valid for physical lattices as well by extending some expansion methods 
developed in the framework of dimer deposition j35| . 

Finally, one would like to extend the boolean CA approach to three-state models 
accounting for the presence of more than one species and, more generally, to models 
displaying complex MF behavior like oscillations |36|-p8[] and phase transitions |39l,f40 
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APPENDIX A: THE SYMMETRIC CPD 

In this symmetric version of the CPD model, a particle at the chosen site looks either 
at the left or at the right neighbor site with equal probability and desorbs if the chosen 
neighbor is occupied. In this case, the dynamical rule reads: 

m(t + At) = m(t) - f$(t) m(t) {&(*) rn-i(t) + (1 - &(t)) n i+1 (t)} , i = l,..., N. 

(Al) 

The additional decision variable is equal to one if the left neighbour is chosen and 
zero otherwise. Clearly, the additional choice between left and right may change the 
global coverage in a realization characterized by a given path {£jv(0), . . . , £jv(*)}- In 
contrast, one easily checks that averaging Eq. (|A1|) leads to the balance equations (|20f) 
and (|26|) valid for the asymmetric system. Clearly, the contributions to the destruction of 
a cluster by reactive events between two particles inside the cluster are the same in both 
models. In the asymmetric model, a cluster may also be destroyed by the choice of its 
rightmost particle, which will react with its right occupied neighbour with probability 
one. In contrast, this event will only take place with probability | in the symmetric 
model, but there is an additional such contribution due to the reaction of the leftmost 
particle in the cluster with its left occupied neighbor. Therefore, both models will lead to 
the same balance equations for the cluster probabilities Pk(t) and Pk(t). This conclusion 
is supported by comparison of MC simulations performed according to the rules (H) and 
( |A1|) . A similar argument can be applied to the symmetric CTD. 



APPENDIX B: DYNAMICS OF PARTICLE ISLANDS 



The definition of a particle cluster introduced in 11 (J was nonexclusive, meaning that 
a cluster of a given size could contain smaller clusters. In the following, we will consider 
a more restrictive definition in which only isolated strings of particles are regarded as 
distinct clusters. These "particle islands" are characterized by the nonvanishing product 



(l-ni_i(t)) 



i+Jfc-l 

n 

j=i 



(1 - n i+k (t)), 



(Bl) 



where k is the size of the island. In a iV-site periodic lattice, the total number of islands 
of size k is given by 



N 



iV/ >fc (*) = E( 1 - n i-i(*)) 



i=l 



i+k-l 

n 



(1 -n i+fc (i)), 



(B2) 



and the total number of islands is Nj(t) = Y^ k =i Ni lk {t) (the largest island can at most 
have N — 1 sites in the ring). With this definition, islands can only be created, never 
destroyed, even though its size may decrease in time. In order to create a new island at 
time t, an existing one must be split up by removal of an internal particle (not at the 
edge of an island). Therefore, the time evolution of Nj(t) will be given by 



N 



Nj(t + At) = N T (t) + 7ii-i(t) ni{t) n l+1 (t) 



(B3) 



i=l 



In the limit t — > oo, the iteration of this formula yields the total number of particles 
Na(oo) in the steady state: 



oo N 



N A (oo) = Nj(oo) = N I (0) + J2T,^N > ( t ) rii(t)n i+1 (t). 



(B4) 



t=o i=l 



The validity of this formula for the asymptotic coverage induced by a single realization 
can be easily checked numerically. 



Following paragraph [II D| , we now consider an ensemble of realizations starting from 
the same initial condition. The probability to find an island of k consecutive sites can 
be defined in a similar way to Pk(t): 



N 



N 



E4 l) W, k = i,...,N-i. 



N 



(B5) 



i=l 



with 



4 i) (t) = (l-^ 1 (t)) n nj(t)(l-n i+k (t)). 



i+k—l 



After some algebra, one obtains 

h(t) = P k {t) - 2P k+1 {t) + P k+2 {t). 



(B6) 



(B7) 



Thus, once the solution for the set of P k {t) is known, it is easy to determine the time 
evolution for the islands. A similar relation is found for coarse-grained initial conditions. 
In the continuous-time limit, one has: 

i+k-l 

h(t)= ( (l - m-tit)) J] n i(t) ( x - m+k(t)) ) 

j=i 

= P k (t) - 2P fc+1 (t) + P k+2 {t) = [1 - 2p exp(-t) + p 2 exp(-2t)] P k (t). (B8) 

Note that I\(t) increases monotonically in time, since the number of isolated particles 
always increases. The situation is less obvious for multiparticle islands (Fig. |9|). If the 
initial coverage is high enough, their number will first be increased due to the breaking of 
larger islands, however they will sooner or later be themselves reduced to single-particle 
islands by the ongoing reactions. 

The situation is slightly different in the CTD, since two-particle islands can indeed 
be destroyed. The evolution of Nj(t) in a single realization is given by: 

N 

n i+ i(t) n i+2 (t) 

8=1 

TV 

(1 - rn-i(t) ) rn(t) m+i(t) (1 - n i+2 {t) ). (B9) 

i=i 

The first term describes the creation of a new island by reaction of two internal particles, 
while the second term stands for the destruction of a two-particle island. The total 
number of particles in the steady state can be written as 

oo N 

N A (oo) = Nj(0) + ££d?(t)n i (t)n m (t) (n^t) + n i+2 {t) - 1). (BIO) 

t=0 i=l 



The dynamics of the islands can be easily determined by using the relations (|B7p and 
(|B8|), which also hold in this case. 
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FIG. 1. Reaction step for the CPD. If the chosen site and a randomly chosen neighbor are 
filled, the reaction proceeds with rate &r. 

FIG. 2. (a) Mean coverage c(t) = P\(t) for the four different initial configurations CI — C4 
computed from MC simulations over 10 4 realizations. Both the dynamics and the asymptotic 
value of the coverage are in excellent agreement with formulae (^) and (|24|). The final satu- 
ration at different coverage values reflects the dependence on the initial conditions, and hence, 
the weak ergodic properties of the system, (b) Time evolution of the coarse-grained coverage 
c(t) = Pi(t) computed from MC simulations over 10 2 realizations for a periodic lattice of 10 4 
sites and p = 0.35,0.5 and 1 (full lattice). As in Fig. |2](a), the final saturation value depends 
on the initial coverage p. Note, in addition, the monotonic dependence for all times of c{t) as 
a function of the initial coverage p. 

FIG. 3. Time evolution of the mean concentration c(t) according to the rate equation (0) 
with kji = 1 (curve MF), the exact solution (|28a| ) (curve ES) and the numerical solution for the 
pair approximation of the truncation hierarchy (^8|) (curve PA) . Initially, the lattice is assumed 
to be fully occupied (p = 1). In contrast to the exact solution, the truncated one decays slowly 
to the MF steady state. The inset shows a log-log plot of the concentration for long times for 
the MF and the truncated solution. 

FIG. 4. Time evolution of /i(i),/2(i) and h(t) computed from MC-simulations over 10 2 
realizations for a periodic lattice with N = 10 4 and p = 0.35. Asymptotically, f%(t) is nearly an 
order of magnitude smaller than fi(t), in accordance with the strong spatial decay predicted 
by formula (|54|). In contrast, the three-point correlation h(t) is of the same order of magnitude 
as f 2 (t). 

FIG. 5. p-dependence of the asymptotic ratio h(cx)/ f 2(00) computed from the formulae 
(|54| ) and (|5q). For decreasing p, three-point fluctuations become increasingly important. 

FIG. 6. Analytical solution of the first four evolution equations for intervals of vacant sites 
(p = 1). The quantities S1-S4 grow monotonically due to the empty segments created by the 
reaction. 

FIG. 7. Reaction step for the CTD. In contrast to the CPD case, both interacting particles 
desorb when the event takes place. 

FIG. 8. (a) Mean coverage c(t) = Pi(t) for a chain ring of 10 2 sites for the four different 
initial configurations CI — C4 computed from MC simulations over 10 4 realizations (c.f. Fig. 
|(a) ). Notice, again, the dependence of the asymptotic coverage on the initial conditions, (b) 
Time evolution of the coarse-grained coverage c(t) = P\(t) computed from MC simulations 
over 10 2 realizations for a periodic lattice of 10 4 sites with p = 0.35,0.5 and 1 (full lattice). 
In contrast to the CPD case, the behavior of c(t) is no longer monotonic for sufficiently large 
times. 



FIG. 9. Time dependence of Ii(t), ^(t) and /3(f) according to formula (|B8| ) for an infinite, 
initially full chain (p=l). 



